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ABSTRACT 

Aims. We study the acceleration and radiation of charged particles in the shock waves of supernova remnants using a recent version of the 
"box model". According to this, particles are accelerated in an energy-dependent region around the shock by the first order Fermi mechanism 
and lose energy through radiation. 

Methods. The particle distribution function is obtained from a spatially averaged kinetic equation that treats the energy losses self-consistently. 
There exists also a second population that consists of those particles that escape behind the shock where they also radiate. The energy 
distribution of this population is calculated in a similar manner. 

Results. The application of the model to the supernova remnant RX J1713. 7-3946, which was recently confirmed as a TeV source by H.E.S.S., 
shows that the X-ray emission can be attributed to electron synchrotron radiation while in y-rays there are contributions from both electrons 
and protons, with protons playing the dominant role. Additionally, there are strong indications that particles diffuse in turbulence that has a 
Kolmogorov spectrum. 
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1. Introduction 

Supernova remnants (SNRs) have long thought to be the sites 
of galactic cosmic ray acceleration, at least up to energies 
close to that of the "knee" (~ 10 15 eV) of their spectrum 
(Krymskii 1977; Bell 1978; Blandford & Ostriker 1978). Direct 
evidence for this came only during the last decade with the 
observation of non-thermal X-ray emission from SNRs. This 
can be best explained as synchrotron radiation by electrons that 
have been accelerated up to energies ~ 100 TeV (Reynolds & 
Keohane 1999). Despite this success, evidence for the presence 
of a relativistic hadronic component is less certain. Accelerated 
hadrons in SNRs are expected to emit high energy gamma rays 
from the decay of neutral pions (tt°) produced in proton-proton 
collisions (Dorfi 1991). These are best observed in the TeV 
regime (Drury et al. 1994), thus SNRs have been one of the 
main targets for the ground based Cerenkov detectors (for a re- 
view see Weekes et al. 2002). However, the above picture is 
complicated from the fact that the same high energy electrons 
which produce the synchrotron X-rays will Compton up-scatter 
the microwave photons to TeV energies (Mastichiadis 1996; 
Pohl 1996), contributing an extra component in this regime. 
Therefore one expects that any TeV detection of a SNR can- 
not be unequivocally attributed to the presence of a relativistic 
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hadronic (or leptonic) component. Instead careful modeling is 
required to determine the origin of these gamma-rays. 

In modeling the multi-wavelength spectrum of SNRs most 
authors assume for the particle energy distribution (electrons or 
protons) the form of a power-law with an exponential cut-off. 
The power-law is a natural consequence of Fermi acceleration 
while the exponential turnover is supposed to arise from energy 
losses or SNR age limitations. Despite the theoretical support 
for these assumptions, this approach is still phenomenological. 
On the other hand, the use of extensive numerical codes has 
proved a valuable tool for studying the acceleration, especially 
as they allow more complex phenomena, such as shock modi- 
fication, to be taken into account (Ellison et al. 2001; Berezhko 
et al. 2002). The application of an analytical model however is 
interesting as it offers a good opportunity to get an insight at 
some of the basic principles of particle acceleration. 

A development in the direction of an analytical approach to 
the problem of particle acceleration has been made recently by 
Drury et al. (1999). This approach has expanded on the concept 
of the so-called "box" model for particle acceleration (for a re- 
view see Kirk et al. 1994) by assuming an energy-dependent 
size for the acceleration region. This inserts some modifica- 
tions on the particle spectra and is thus interesting to apply this 
model to fit the SNR multi-wavelength emission. A good start- 
ing point for this is the first SNR with confirmed TeV emission, 
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RX J1713.7-3946 (Muraishi et al. 2000; Enomoto et al. 2002; 
Aharonian et al. 2004, 2006). 

RX J1713.7-3946 (G347.3-0.5) is a shell-type supernova 
remnant that was discovered in X-rays during the ROSAT All 
Sky Survey (Pfeffermann & Aschenbach 1996). It is a young 
SNR with an age of ~ 2000 years, possibly associated with the 
supernova that took place in AD393 (Wang et al. 1997) and is 
located in the Galactic plane at a distance ~ 1 kpc. The X-ray 
emission of RX J1713.7-3946 is purely non-thermal in nature 
(Koyama et al. 1997; Slane et al. 1999) and its X-ray morphol- 
ogy, as revealed by Chandra, is quite complicated consisting 
of bright filaments and dark voids (Uchiyama et al. 2003) in- 
dicating spatially inhomogeneous magnetic field and/or den- 
sity of relativistic particles. In y-rays, the remnant was origi- 
nally observed by the CANGAROO telescope (Muraishi et al. 
2000; Enomoto et al. 2002) and recently confirmed as a high 
energy source by H.E.S.S. (Aharonian et al. 2004, 2006). The 
remnant also exhibits a weak radio emission (Lazendic et al. 
2004). The broadband energy spectrum of RX J1713.7-3946 
has been modeled by many authors, assuming either a lep- 
tonic (Uchiyama et al. 2003; Pannuti et al. 2003; Lazendic 
et al. 2004) and/or a hadronic origin (Enomoto et al. 2002; 
Aharonian et al. 2006; Berezhko & Volk 2006) for the y-ray 
emission, but without a conclusive result so far. 

In this paper we apply the Drury et al. (1999) model to the 
supernova remnant RX J1713. 7-3946 and draw some conclu- 
sions about its properties. In section 2 we review the main char- 
acteristics of the model and obtain the distribution function of 
accelerated particles as well as that of the particles which have 
escaped behind the shock. In section 3 we derive the expected 
multi-wavelength photon spectrum from these distributions and 
apply it to RX J1713. 7-3946. Finally, in section 4 we summa- 
rize and discuss our results. 

2. The acceleration model 

Charged particles interacting with the shock wave of a su- 
pernova remnant can be accelerated to high energies through 
first order Fermi acceleration. Alfven waves in the background 
medium scatter the charged particles, changing their direction 
but not their energy. Particles are thus confined to a region 
around the shock and are forced to cross it many times. The dif- 
ference in the gas velocity upstream and downstream results in 
the statistical acceleration of particles after many shock cross- 
ings. While confined around the shock particles also lose en- 
ergy because they emit radiation. At the same time a fraction of 
them escapes behind the shock towards the interior of the rem- 
nant where it also radiates. We next calculate the distribution 
function of these two populations when the accelerated parti- 
cles are electrons as well as when they are protons. We denote 
quantities in the acceleration region with the subscript 1 and in 
the escape region with the subscript 2. 

Although the shock in a SNR is spherical we can treat the 
acceleration in the plane shock geometry as long as the size of 
the acceleration region is much smaller than the SNR radius. 
We assume the presence of a uniform magnetic field that is 
normal to the shock. We also adopt the test particle approxima- 
tion, i.e. we assume zero particle pressure. With these assump- 



tions, acceleration takes place in a region around the shock of 
size L = L u + L d = K u /U u + K d /U d where K u(A) , U u(d) are 
the upstream (downstream) diffusion coefficient and gas veloc- 
ity respectively. Assuming that the diffusion coefficient is con- 
tinuous across the shock, i.e. K u = — K, the size of the 
acceleration region can be further written L = K(r + 1)/U U 
where r = U u /Ud is the compression ratio of the shock. 
Following Drury et al. (1999) we let the diffusion coefficient 
to be momentum-dependent and adopt a general power-law 
dependence K(p) = Kp s with 6 > 0. This choice makes the 
mathematical analysis easier but also covers some physically 
interesting cases. For example, the case 6 = 1 corresponds 
to Bohm diffusion, 6 = 1/3 to diffusion from Alfven waves 
with a Kolmogorov spectrum of turbulence and 8 = 1/2 to a 
Kraichnan spectrum. The same power-law dependence holds 
for the size of the acceleration region L(p) = K(r + l)p s /U u as 
well as the acceleration timescale which is given by the relation 
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If the accelerated particles are electrons they also experi- 
ence energy losses through synchrotron radiation and/or in- 
verse Compton scattering (IC). In the case where the latter is 
in the Thomson limit, we can write for the particles total losses 
p = -aip 2 , where 



4 cr T (f/ Bl + t/ph) 
a\ = t — 



(2) 



is essentially the sum of the energy densities of the magnetic 
field f/ B , and the target photon field f/ ph , <x T is the Thomson 
cross section, m the electron mass and c the speed of light. If 
there were no energy losses, particles would escape from the 
acceleration region with the downstream flow at the bulk veloc- 
ity Ud- The presence of losses however generates a downward 
flux in momentum space and when combined with the increas- 
ing size of the acceleration region gives an additional escape 
process. The velocity of escape thus increases and the escape 
timescale is given by 
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The second term in the denominator is the addition of the Drury 
et al. (1999) model to the standard box model and vanishes for 
constant diffusion coefficient or no energy losses. 

The evolution of the spatially averaged isotropic distribu- 
tion function of particles at the shock front, N\{p,f), is gov- 
erned by the continuity equation 



dN x d 
~~dt~ + dp 
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The first term in brackets describes acceleration at the momen- 
tum dependent rate r^ c and the second the synchrotron-type 
energy losses. The last term in the left hand side describes par- 
ticle escape at the momentum dependent rate and the source 
term Q(p, t) = Q 8(p - p )fl-(t) represents injection of particles 
at low momenta p a with the constant rate Q ; 5{x) and i?(x) are 
the Dirac delta and Heavyside (or unit-step) functions respec- 
tively. 



K. Moraitis and A. Mastichiadis: A two-zone model for the emission from RX J1713.7-3946 



3 



From the equation describing the evolution of a particle's 
momentum 



dp 
dt 



(5) 



it follows that particles can be accelerated up to the momen- 
tum where the acceleration and cooling timescales are equal, 
namely 
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Particles reach this maximum allowed momentum if given 
enough time, otherwise their momentum is determined from 
the solution of Eq. (5). This equation can be solved analyti- 
cally for the time a particle of initial momentum p Q needs to 
reach momentum p and the solution is 



hiP) = fcool.l 
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where jF\ is the hypergeometric function and ? coo i_i = 
(fi/'max) 1 is the minimum cooling timescale. The inverse 
function of t\(p) gives the maximum momentum of a particle 
at time t and is denoted as p\it). This maximum momentum in- 
creases with time until energy losses become important and the 
acceleration saturates at p max as can be seen in Fig. 1. The time 
required for this (in units of ? C ooi,i) is a decreasing function of 
the index 5. 

The solution of Eq. (4) for the distribution function is then 
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in the momentum range p Q < p < p\(f) and zero otherwise; f a>0 
is the minimum acceleration timescale, t. dcc (p ), and the index 
b = + j^j. The accelerated particles form a power- 

law of index si = - 8, smaller by 6 from the index pre- 
dicted by the standard box model, that extends from the in- 
jection momentum p up to the maximum momentum p\if). 
Given enough time the particle distribution function reaches 
steady-state and is then given by Eq. (8) in the momentum 
range p < p < p max . In steady-state the power-law gets mod- 
ified since the first term in brackets in Eq. (8) is then zero at 
P = Pmax and a pile-up or a cut-off of particles forms depend- 
ing on whether the index b is positive or negative respectively. 

Particles that escape towards the interior of the remnant 
play the most important role in the emission (Ball & Kirk 1992; 
Aharonian & Atoyan 1999) and we calculate next their en- 
ergy distribution. The distribution function of escaping parti- 
cles obeys a similar equation to Eq. (4) where now the source 
term is the escape rate of accelerated particles while there are 
no acceleration or escape terms 



IT - T P ( aipN2 ) 



Ndp,t) 

fesc(p) ' 



(9) 




Fig. 1. The evolution of the upper limit of the particle distribution 
function p\{t) (solid line) for parameters p a = 0, 6 = 0.3 and 
B t = 10 pG. Time is measured in units of the cooling timescale 
^cooi.i = 2200 (p m ax c/50 TeV)~'yr. The evolution of the cooling mo- 
mentum p coo i(t) for the cases B2 = 10 pG (dashed line), B2 = 5 pG 
(dotted line) and B 2 = 2 pG (dash-dotted line) is also shown. 



The energy losses have the same form as in the acceleration 
zone (see Eq. 2) but depend on the magnetic field of the es- 
cape zone, Z?2- The solution of Eq. (9) can be easily obtained 
numerically (and semi-analytically) but we prefer a simplified 
analytic approach. If the energy losses are negligible the solu- 
tion is simply the escape rate integrated over time 
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which, for p <K piit), is a power-law with the standard box 
model index S2 = si + 6 = ^y. This result remains valid 
when losses are present but only for those particles that had 
not enough time to cool, i.e. those with momentum less than 
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This equation describes the evolution of momentum of the 
highest energy particles as they cool in the escape zone and 
is shown in Fig. 1. On the other hand, particles with momen- 
tum exceeding p coo \it) had the necessary time to cool and their 
distribution is steeper from the injected one by one power of p, 
i.e. 



N 2 ,cooiip, * N 2 ip, t) 



Pcoolit) 



(12) 



The escaping particles thus have a broken power-law distribu- 
tion, given by Eq. (10) for p < p coo \it) and by Eq. (12) for 
Pcoolit) < p < piit). The break occurs at the cooling mo- 
mentum which depends on the magnetic field of the escape 
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region. The extent of the cooled part of the particle distribu- 
tion, pi(t) - p C ooi(t), is greater for higher values of Bi as is 
shown in Fig. 1. The steady-state escape distribution consists 
only of cooled particles and hence has a power-law index equal 

to S2 + 1 . 

We next consider the case that the accelerated particles are 
protons. Since protons lose substantial energy only through the 
production of pions and the corresponding cooling timescale 
is much larger than the age of the remnant, f pp 5 10 15 s » 
?sn - 2000 yr, we can neglect their energy losses. We thus 
solve Eqs. (4) and (9) without the energy loss term and for a 
source function that provides protons of momentum p Q at the 
rate Q a . The proton energy distribution in the acceleration zone 
has then the simple form 



Ni(p,t) 



Qok,o I P 



Po \Po 
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in the momentum range p < p < p\(f) and f a ,o = td CC (p ). The 
upper limit of the proton distribution follows from the solution 
of Eq. (5) with a\ = and is given by 



(14) 
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The distribution in the escape region is 
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in the momentum range p a < p < p\(t), where again t\(p) is the 
inverse function of pi(t). We note that the power-law indices si, 
S2 are the same as in the leptonic case. 

3. Application to the supernova remnant 
RX J1 71 3.7-3946 

We now apply the model developed in the previous section to 
the supernova remnant RX J1713. 7-3946. We first consider 
the case where only electrons are responsible for the observed 
emission, that is we neglect the possible contribution from a 
hadronic component. As discussed in section 2 there exist two 
relativistic electron populations, one in the acceleration region 
and one in the escape region behind the shock. These elec- 
trons will produce synchrotron radiation in the magnetic field 
of the two regions, assumed equal for simplicity. Electrons will 
also Compton up-scatter the photons of the cosmic microwave 
and infrared background radiation fields which are the same in 
both zones. For the microwave photon field we take the stan- 
dard values of temperature T M w = 2.7 K and energy density 
f/MW = 0.25 eV crrT 3 while for the infrared field we assume 
T IR = 25 K and Ujr = 0.05 eV cnr 3 (Aharonian et al. 2006). 
In the possible emission processes we also include non-thermal 
bremsstrahlung radiation and assume for the ambient hydrogen 
density the average value n ~ 1 crrT 3 . 

The photon spectrum derives by folding the two electron 
distribution functions with the single particle emissivity for the 
synchrotron, IC scattering and bremsstrahlung processes and 
thus depends on the same parameters that determine the elec- 
tron distribution function in the acceleration and escape zones. 



These parameters are the compression ratio of the shock r, the 
index of the diffusion coefficient 5, the magnetic field B, the 
maximum possible electron energy £ e ,max = PmaxC and the nor- 
malization constant Q Q and can be determined from the fit of 
the model to the observed spectrum. The injection momentum 
of electrons does not alter the spectrum significantly and so we 
take the arbitrary value p D = 5 MeV/c. We should also note 
that the calculation of the photon spectrum is for the present 
day t = ?sn- 

The multi-wavelength spectrum of RX J1713. 7-3946 is 
shown in Fig. 2. Radio data from the northwest part of the 
remnant taken with ATCA (Lazendic et al. 2004) are shown 
here multiplied by a factor of 2 to cover the whole remnant. 
Likewise, the X-ray data of ASCA (Enomoto et al. 2002) are 
multiplied by a factor of 4 so that the total X-ray flux from 
the remnant is obtained (Tomida 1999). Finally, the TeV data 
are the recent observations of H.E.S.S. for the entire rem- 
nant (Aharonian et al. 2006). We also plot the EGRET data 
of the nearby source 3EG J1714-3857, since as emphasized 
by Reimer & Pohl (2002) they set an upper limit to the GeV 
emission of RX J1713.7-3946. 

We first examine a purely leptonic model in which the 
X-ray emission is due to synchrotron radiation and the y-ray 
emission is due to IC scattering. From the observational data 
plotted in Fig. 2 we deduce that the X-ray emission peaks at 



1 keV and the y-ray emission at s 



.TeV 
max 



- 1 TeV. Since 

the peak frequencies scale as El m . dx B and i^max respectively 
(IC scattering takes place mostly in the Thomson regime), we 
find that the magnetic field should be B ~ 40 pG and also that 
£e,max ~ 45 TeV. For these values, electron cooling is impor- 
tant and the escaping particles have a broken power-law dis- 
tribution which gives the correct shape for the extended y-ray 
spectrum. However, the X-ray flux level is then overestimated 
by a factor of 15 and one should use the so-called magnetic 
filling factor to get the correct flux level. This however can- 
not be done in our model since the emitted radiation is cal- 
culated self-consistently from the particle losses. We thus see 
that a synchrotron-IC model cannot explain at the same time 
both the relative flux levels and the relative position of the peak 
frequencies. The possibility that the y-ray emission is due to 
bremsstrahlung radiation can most probably be ruled out, since 
the ambient hydrogen density required to fit the spectrum is 
n ~ 300 crrT 3 , a value that might be too high. We conclude that 
the observed radiation from this remnant cannot be of purely 
leptonic origin. 

We thus consider that protons are also accelerated in 
the remnant. As a secondary electron population produced 
from charged pion decay cannot explain the X-ray emission 
(Mastichiadis 1996), primary electrons are again required. 
However, in this case both electrons and protons contribute to 
the y-ray regime. The inclusion of a relativistic proton popu- 
lation introduces one more free parameter, the energy content 
in protons S p (related to Q ). For the minimum proton energy 
we adopt the value £p >mm = 2 GeV while the maximum pro- 
ton energy is determined from Eq. (14). For the emission from 
neutral pions we use a 5-function approximation. 

We fit our model to the multi-wavelength spectrum of 
RX J1713. 7-3946 for the parameters r, 6, B, £e,max, Qo and 



K. Moraitis and A. Mastichiadis: A two-zone model for the emission from RX J1713. 7-3946 



5 



<2o- The resulting theoretical spectrum is shown in Fig. 2 to- 
gether with the observational data. Although there is a small 
disagreement with the first two H.E.S.S. points (at 215 GeV 
and 255 GeV), the fit is quite satisfactory and has a chi square 
X 1 — 23 with 22 degrees of freedom (when the last H.E.S.S. 
point at 35 TeV is excluded from the fit). From the derived 
values of r = 3.8 and 8 = 0.3 we can calculate the two par- 
ticle power-law indices s\ = 1. 11 and S2 = 2.07. We note 
that the value of r implies a strong shock and the value of 6 
a nearly Kolmogorov-type diffusion K(p) oc p 03 . The mag- 
netic field has the value B = 15 yuG. From this value and 
the maximum electron energy Ee.max = 100 TeV we find that 
the cooling timescale is f coo i = 500 yr and hence the en- 
ergy losses in the escape region are marginally important. We 
also deduce that electrons in the acceleration zone are close 
to steady-state, /?i(/sn) = 0.7 p m - ix , while those in the escape 
zone have the broken power-law distribution with the break at 
Pcooi(^sn) = 20 TeV/c = 0.2 p max . For the protons the maxi- 
mum energy at the present day is 160 TeV. The normalization 
constant Q Q = 6 10 40 s _1 determines the total energy content 
in electrons from both zones £ e = 3 10 47 erg. The proton nor- 
malization constant is Q Q = 2 10 41 s _1 and the correspond- 
ing energy content in protons is £ p = 2 10 50 erg, dominating 
the total energy content in relativistic particles, i.e. £ tot - fi p . 
The majority of this energy is carried by particles in the es- 
cape region. For the typical value of supernova explosion en- 
ergy of a few times 10 51 erg we find that a fraction ~ 10% 
of it is needed in order to explain the non-thermal emission of 
RXJ1713.7-3946. 

From these parameters we find that the diffusion coefficient 
has the form 

K {p) = 3 10 26 ( ^ r) 2 ( P - X 3 cm 2 s- 1 

yiJ l610 3 kms-V \ 100 TeV c- 1 / 

and is greater than the lower limit of Bohm diffusion for all 
particle momenta less than 

If we thus adopt for the shock velocity the value U s h = 
6 10 3 km s _1 then the diffusion coefficient we derive does not 
violate the lower limit of Bomh diffusion for any value of par- 
ticle momentum in our model. 

4. Summary and Discussion 

In the present paper we have attempted a fit to the 
multi-wavelength spectrum of the supernova remnant 
RX J1713.7-3946 based on the box model of Drury et al. 
(1999). According to this, particles diffuse in an energy- 
dependent region around the shock front of the remnant and 
are accelerated by the first order Fermi mechanism. The 
particle distribution functions (both for electrons and protons) 
come as solutions to two coupled partial differential equations, 
one characterizing the particles inside the acceleration region 
and the other the ones that have escaped from it. Consequently, 
one can calculate the radiated photon spectra by folding 
the derived particle distribution functions with the single 
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Fig. 2. The multi-wavelength spectrum of RX J1713.7-3946 together 
with the best-fit model with parameters given in text. We also show 
the GeV upper limits from the nearby source 3EG J1714-3857. The 
dashed line represents the emission from neutral pions, the dotted line 
from IC scattering and the solid line is the total emission. Also plotted 
is the bremsstrahlung emission for the average ambient density n = 
1 cirr 3 as the dash-dotted line. 



particle emissivity for the radiative processes assumed. For 
the case of electrons the corresponding radiative processes are 
synchrotron, IC scattering and bremsstrahlung emission, while 
protons produce radiation only through the decay of pions. 

We found that, within the strict assumptions of the box 
model, a purely leptonic origin of the radiation cannot explain 
the observations. The leptonic model is viable only in the case 
where the magnetic field in the escape region has a value large 
enough as to allow for substantial particle cooling during the 
lifetime of the remnant. The escaping particles then have a bro- 
ken power-law distribution that can explain the extended y-ray 
spectrum of RX J1713.7-3946. For such a high value of the 
magnetic field though, the level of X-ray emission exceeds sig- 
nificantly the observed one. We thus reach to the same con- 
clusion with all the other leptonic models used for this source, 
i.e. that we cannot explain at the same time the relative posi- 
tion of the synchrotron and IC peaks and their relative fluxes, 
unless we use the filling factor. We should emphasize however 
that the above conclusions are valid under the assumptions of 
a spatially averaged model and a plane-parallel shock geome- 
try. In a model where the spatial dependence is treated in de- 
tail one could include other effects like variable injection rate 
at different supernova radii and/or expansion losses, although 
the error from the latter should not be important. Whether a 
purely leptonic model would be more successful in such a case 
is something that remains to be investigated. 

The inclusion of the relativistic hadronic population pro- 
vides an additional component in y-rays and a satisfactory fit 
can then be obtained. However, the hadronic component can 
be treated only approximately in our model since the test par- 
ticle approximation that we used is not strictly valid and the 
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shock gets modified by the hadrons. To take this into account 
one should use calculations based on numerical models. 

From the fit to the spectrum of RX J1713.7-3946 we also 
deduced the dependence of the diffusion coefficient on energy. 
We found that a nearly Kolmogorov-type diffusion coefficient 
(6 = 0.3) can reproduce the observed spectrum quite well. The 
lcr estimate for the parameter 6 is well concentrated around 
this value, 0.1 < 6 < 0.5. This is because a change in 6 mostly 
influences the maximum electron and proton energies and these 
are well defined from the observed spectrum. We note that this 
range of 6 allows marginally for a Kraichnan-type turbulence 
spectrum, while the case of Bohm diffusion is not favored by 
our results. 

We conclude that the model we used, despite certain as- 
sumptions, has two main advantages over other models, in that 
a) the particle distribution functions are calculated from a the- 
oretical model and are not assumed a priori and b) it can give 
a method of estimating the energy dependence of the diffusion 
coefficient in supernova remnants. 
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